# multinet.mod    
# Yu Liu, 01/2006, model multi layer network all together.


# layout.mod
# by Yu Liu, 03/2001.   Based on arcflow.mod
# Arc-flow model for topology and working capacity allocation 
# considering all single NODE failures

# Declare the bottom layer topology with VERTS and EDGES
set VERTS;
set EDGES within {i1 in VERTS, VERTS diff {i1}};

# Flows on the bottom layer is FLOWSL where suffix 'L' for "Lower Layer"
set FLOWSL := {i1 in VERTS, VERTS diff {i1}};
# Each flow has a bandwidth requirement in parameter ML, default at 0 for bottom layer
param ML {FLOWSL} default 0 >=0; 

# matrix oBL has element 1 when the vertex is the original of the edge, 0 otherwise. 
param oBL {VERTS,EDGES} binary default 0;
# matrix dBL has element 1 when the vertex is the destination of the edge, 0 otherwise.
param dBL {VERTS,EDGES} binary default 0;

# matrix BL is the vertex to edge incidence matrix at bottem layer
# The element is 1 when the vertex is the original of the edge, 
#               -1 when the vertex is the destination of the edge,
#                0 otherwise.
param BL {VERTS,EDGES} integer default 0 >= -1 <= 1;

# matrices oDL, dDL, and DL are the flows to vertex incidence matrices at bottom layer
# They are similar to the declarition of oBL, dBL, and BL
param oDL {FLOWSL,VERTS} binary default 0;
param dDL {FLOWSL,VERTS} binary default 0;
param DL {FLOWSL,VERTS} integer default 0 >= -1 <= 1;


# Declare top layer topology
# node set is a subset of lower layer vertices, links are 
set NODES within VERTS;
set LINKS within {i1 in NODES, NODES diff {i1}};

# Declare flows on top layer
set FLOWS := {i1 in NODES, NODES diff {i1}};

# Each flow has a bandwidth requirement in parameter M
param M {FLOWS} default 1 >=0; 

param oB {NODES,LINKS} binary default 0;
param dB {NODES,LINKS} binary default 0;
param B {NODES,LINKS} integer default 0 >= -1 <= 1;
param oD {FLOWS,NODES} binary default 0;
param dD {FLOWS,NODES} binary default 0;
param D {FLOWS,NODES} integer default 0 >= -1 <= 1;

# matrix H maps each top layer LINK into a path on the bottom layer EDGES
var H {LINKS, EDGES} binary default 0;

# FIXME two auxilary matrices used for what? They are declasred as variables...
#var T1 {LINKS} binary default 1;
#var T2 {(i1,i2) in LINKS, NODES diff {i1,i2}} binary default 1;
#var A {FLOWS, LINKS} binary;
#var w {LINKS} integer default 0;

# total numbers of nodes and links on both layers
param MaxVert integer;
param MaxEdge integer;
param MaxNode integer ;
param MaxLink integer ;

# build a set of CUTS, 
# each of its elements hold a set of nodes.
# 
set POWNODES := 1 .. 2**MaxNode - 2; #get rid of 0 and 2^N -1, null cuts
set CUTS {k in POWNODES} := {i in NODES: (k div 2**(i-1)) mod 2 = 1};

# optimize the cost of top layer links on bottom layers, 
# assuming the unit link cost on the bottom layer
minimize topo_cost: 
  sum {(i1,i2) in LINKS, (j1,j2) in EDGES} H[i1,i2,j1,j2];

# Interlayer mapping matrix maps top layer links as the paths at bottom layer.
# These paths fulfill the flow conservation (mass balance) constraints
# The following formula assumes the first group of vertices (at bottom layer) are nodes (at top layer).
# the rest vertices are only for the bottom layer, so their net in/out flows are zeroes.

# H {B^l}^T = [ B^T | 0 ] 
s.t. mass_baH {(i1,i2) in LINKS, n1 in NODES}:
  sum {(j1,j2) in EDGES} BL[n1,j1,j2] * H[i1,i2,j1,j2]
      = B[n1,i1,i2];
s.t. mass_baH2 {(i1,i2) in LINKS, n1 in VERTS diff NODES}:
  sum {(j1,j2) in EDGES} BL[n1,j1,j2] * H[i1,i2,j1,j2]
      = 0;

# yliu: 01/28/06: symH gives smaller solution space, that is why not used earlier
# used it now for published results
# added by Korn, H should have the same mapping for both directions
s.t. symH {(i1,i2) in LINKS, (j1,j2) in EDGES}:
	H[i2,i1,j2,j1] = H[i1,i2,j1,j2];


# This formula assures the single failure of bottom layer will not partition the top layer.
# It is the matrix format of the 

# CH < Ce  or  C(e-H)>0
s.t. survH {(j1,j2) in EDGES, k in POWNODES}:
  sum{v1 in CUTS[k], v2 in NODES diff CUTS[k]:(v1,v2) in LINKS}
    (1- H[v1,v2,j1,j2]) >= 1; 

# problem to find interlayer mapping matrix H 
# that assures that the top layer is resilient to single link failure at the bottom layer
problem designH: topo_cost, mass_baH, mass_baH2, survH, symH,  H;

# Where those parameters are used? for what?  yliu:01/2006
param MaxFlow integer default 10;
param MaxHop integer default 60;

# dummy finitions for reading work/backup path
param t1 integer;
param t2 integer;

# workcap.mod
#

param c{LINKS} default 1;

var P {FLOWS, LINKS} binary default 0;
var w {LINKS} default 0;

# min c^t w
minimize w_cost:  sum {(i1,i2) in LINKS} c[i1,i2]*w[i1,i2];

# P B^T = D
s.t. mass_baW {(r1,r2) in FLOWS, n1 in NODES}:
  sum {(i1,i2) in LINKS} P[r1,r2,i1,i2]*B[n1,i1,i2] 
      = D[r1,r2,n1];

s.t. cap_aggrW {(i1,i2) in LINKS}:
  w[i1,i2] = sum{(r1,r2) in FLOWS} P[r1,r2,i1,i2]*M[r1,r2];

s.t. work_sym {(r1,r2) in FLOWS, (i1,i2) in LINKS}:
  P[r1,r2,i1,i2] = P[r2,r1,i2,i1];
# P is symatric for symetric flows

problem find_work:  w_cost,  mass_baW, cap_aggrW, work_sym, P, w;


# sparecap.mod    
# by Yu Liu, 07/04/2001.
# Arc-flow model for spare capacity allocaiton

set FAILS := 1..MaxEdge;
param F {FAILS, LINKS} binary default 0;
# need to be declared somewhere
# for link, node, or arbitrary failure
param U {FLOWS, FAILS} binary default 0;
param T {FLOWS, LINKS} binary default 0;
# can be calculate from F and P

var Q {FLOWS, LINKS} binary default 0;
var G {LINKS, FAILS} default 0;
var s {LINKS} default 0;

minimize s_cost: sum {(i1,i2) in LINKS} c[i1,i2]*s[i1,i2];

s.t. cap_aggrS {(i1,i2) in LINKS, k in FAILS}:
  s[i1,i2] >= G[i1,i2,k];
# s >= G  or s = max G

s.t. spm_comp {(i1,i2) in LINKS, k in FAILS}:
  G[i1,i2,k] = sum{(r1,r2) in FLOWS} 
     M[r1,r2] * (Q[r1,r2,i1,i2] * U[r1,r2,k]);
# G = Q^T M U

s.t. fail_disj {(r1,r2) in FLOWS,(i1,i2) in LINKS}:
  T[r1,r2,i1,i2] + Q[r1,r2,i1,i2] <= 1;
# T + Q <= 1

s.t. mass_baS {(r1,r2) in FLOWS, n1 in NODES}:
  sum {(i1,i2) in LINKS} Q[r1,r2,i1,i2]*B[n1,i1,i2]=D[r1,r2,n1];
# Q B^T = D

s.t. backup_sym {(r1,r2) in FLOWS, (i1,i2) in LINKS:i1<i2}:
  Q[r1,r2,i1,i2] = Q[r2,r1,i2,i1];
# Q is symatric for symetric flows

problem find_spare:
 s_cost, cap_aggrS, spm_comp, fail_disj, mass_baS, backup_sym, Q, G, s;


# sparecap3k.mod
# change from sparecap3.mod is the problem command in the last two lines    
# by Yu Liu, 07/04/2001.
# Arc-flow model for spare capacity allocaiton

# for Model 3 on multilayer networks
# model sparecap.mod;

param cl{EDGES} default 1;

var sl {EDGES} default 0;
minimize sl_cost: sum {(i1,i2) in EDGES} cl[i1, i2]*sl[i1,i2];

s.t. cap_aggrSl {(j1,j2) in EDGES, k in FAILS}:
  sl[j1,j2] >= sum{(i1,i2) in LINKS} H[i1,i2,j1,j2] * G[i1,i2,k];
# sl >= H^T G  or sl = max (H^T G)

problem find_spare3:
 sl_cost, cap_aggrSl, spm_comp, fail_disj, mass_baS, backup_sym, Q, G,
sl;


# workcap_bottom.mod
#
var Pb {FLOWSL, EDGES} binary default 0;
var wl {EDGES} default 0;

# min c^t w
minimize wl_cost:  sum {(i1,i2) in EDGES} cl[i1,i2]*wl[i1,i2];

# P B^T = D
s.t. mass_baW_b {(r1,r2) in FLOWSL, n1 in VERTS}:
  sum {(i1,i2) in EDGES} Pb[r1,r2,i1,i2]*BL[n1,i1,i2] 
      = DL[r1,r2,n1];

s.t. cap_aggrW_b {(i1,i2) in EDGES}:
  wl[i1,i2] = sum{(r1,r2) in FLOWSL} Pb[r1,r2,i1,i2]*ML[r1,r2];

s.t. work_sym_b {(r1,r2) in FLOWSL, (i1,i2) in EDGES}:
  Pb[r1,r2,i1,i2] = Pb[r2,r1,i2,i1];
# P is symatric for symetric flows

problem find_work_b:  wl_cost,  mass_baW_b, cap_aggrW_b, work_sym_b, Pb, wl;

# sparecap_bottom.mod
#
# sparecap_bottom.mod    
# by Yu Liu, 01/22/2006.
# Arc-flow model for spare capacity allocaiton at bottom layer

##already defined in sparecap.mod
##set FAILS := 1..MaxEdge;
##param F {FAILS, LINKS} binary default 0;
# need to be declared somewhere
# for edge, node, or arbitrary failure
param UL {FLOWSL, FAILS} binary default 0;
param TL {FLOWSL, EDGES} binary default 0;
# can be calculate from F and P

var Qb {FLOWS, EDGES} binary default 0;
var Gb {EDGES, FAILS} default 0;

s.t. cap_aggrS_b {(i1,i2) in EDGES, k in FAILS}:
  sl[i1,i2] >= Gb[i1,i2,k];
# s >= G  or s = max G

s.t. spm_comp_b {(i1,i2) in EDGES, k in FAILS}:
  Gb[i1,i2,k] = sum{(r1,r2) in FLOWSL} 
     ML[r1,r2] * (Qb[r1,r2,i1,i2] * UL[r1,r2,k]);
# G = Q^T M U

s.t. fail_disj_b {(r1,r2) in FLOWSL,(i1,i2) in EDGES}:
  TL[r1,r2,i1,i2] + Qb[r1,r2,i1,i2] <= 1;
# T + Q <= 1

s.t. mass_baS_b {(r1,r2) in FLOWSL, n1 in NODES}:
  sum {(i1,i2) in EDGES} Qb[r1,r2,i1,i2]*BL[n1,i1,i2]=DL[r1,r2,n1];
# Q B^T = D

s.t. backup_sym_b {(r1,r2) in FLOWSL, (i1,i2) in EDGES:i1<i2}:
  Qb[r1,r2,i1,i2] = Qb[r2,r1,i2,i1];
# Q is symatric for symetric flows

problem find_spare_b:
 sl_cost, cap_aggrS_b, spm_comp_b, fail_disj_b, mass_baS_b, backup_sym_b, Qb, Gb, sl;

# sparecap2L.mod
# change from sparecap3k.mod  
# by Yu Liu, 01/22/2006.
# Arc-flow model for spare capacity allocaiton for two layer networks
# target for solve model for two layer SCA

# for Model 3 on multilayer networks
# model sparecap.mod;

# include model for bottom layer working
#model workcap_bottom.mod;
#include preworkb.sa
#solve find_work_b;

# include model for bottom layer
#model sparecap_bottom.mod 

s.t. spm_comp_bottom {(j1,j2) in EDGES, k in FAILS}:

  Gb[j1,j2,k] = sum{(r1,r2) in FLOWSL} 

    ( ML[r1,r2] * (Qb[r1,r2,j1,j2] * UL[r1,r2,k]) )

	+ sum{(i1,i2) in LINKS} (H[i1,i2,j1,j2] * sum{(r1,r2) in FLOWS} 

     (M[r1,r2] * (Q[r1,r2,i1,i2] * U[r1,r2,k]) ) ) ;

# Gb = Q^T M U + H ^T (Qb^T Mb Ub)

problem find_spare4:
 sl_cost, cap_aggrS_b, spm_comp, spm_comp_bottom, fail_disj, mass_baS, backup_sym, 
 Q, G, Qb, Gb, sl;


